{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "5e85aaa7-6782-4629-84e0-71a6b7eaa505",
"metadata": {
"nbsphinx": "hidden"
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"%config InlineBackend.figure_format = 'svg'\n",
"\n",
"import warnings \n",
"warnings.filterwarnings(\"ignore\") #ignore some matplotlib warnings\n",
"\n",
"import numpy as np\n",
"\n",
"from mpi4py import MPI\n",
"\n",
"from triqs.plot.mpl_interface import plt, oplot, oplotr, oploti\n",
"plt.rcParams[\"figure.figsize\"] = (6,4) # set default size for all figures"
]
},
{
"cell_type": "markdown",
"id": "e51b0091-8880-4691-9ad1-7cddcff025bd",
"metadata": {},
"source": [
"# Step-by-step guide"
]
},
{
"cell_type": "markdown",
"id": "549bb71a-f72c-47a5-9ba0-1f512be65465",
"metadata": {},
"source": [
"## Step 1 - Impurity many-body Hamiltonian\n",
"\n",
"First we have to construct the local many-body Hamiltonian (`H_loc`) of the Anderson impurity model we want to solve."
]
},
{
"cell_type": "markdown",
"id": "534105f9-15b8-4631-af12-336fcfa9470e",
"metadata": {},
"source": [
"### Examples\n",
"\n",
"- For the **single band Hubbard model** the interaction is density-density only and can be constructed using the Triqs second-quantization operators found in `triqs.operators.*`."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "a1eac349-5da2-4675-9659-417357d3d0f7",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"-1*c_dag('do',0)*c('do',0) + -1*c_dag('up',0)*c('up',0) + 2*c_dag('do',0)*c_dag('up',0)*c('up',0)*c('do',0)\n"
]
}
],
"source": [
"from triqs.operators import n\n",
"\n",
"U = 2.0\n",
"mu = U / 2 # Chemical potential for half-filling\n",
"\n",
"N_tot = n('up', 0) + n('do', 0) # Total density operator\n",
"\n",
"H_loc = U * n('up',0) * n('do', 0) - mu * N_tot\n",
"\n",
"print(H_loc)"
]
},
{
"cell_type": "markdown",
"id": "dedae90a-a720-47cf-b67c-5ae4482bf198",
"metadata": {},
"source": [
"- For multi-orbital models one often use the **Kanamori interaction** which can be built using tools in `triqs.operators.util.*` as follows."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "f55b959b-4469-46d7-9d95-336acdc14e05",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"-7.5*c_dag('do',0)*c('do',0) + -7.5*c_dag('do',1)*c('do',1) + -7.5*c_dag('do',2)*c('do',2) + -7.5*c_dag('up',0)*c('up',0) + -7.5*c_dag('up',1)*c('up',1) + -7.5*c_dag('up',2)*c('up',2) + 2.2*c_dag('do',0)*c_dag('do',1)*c('do',1)*c('do',0) + 2.2*c_dag('do',0)*c_dag('do',2)*c('do',2)*c('do',0) + 0.8*c_dag('do',0)*c_dag('up',0)*c('up',2)*c('do',2) + 0.8*c_dag('do',0)*c_dag('up',0)*c('up',1)*c('do',1) + 4.6*c_dag('do',0)*c_dag('up',0)*c('up',0)*c('do',0) + 3*c_dag('do',0)*c_dag('up',1)*c('up',1)*c('do',0) + 0.8*c_dag('do',0)*c_dag('up',1)*c('up',0)*c('do',1) + 3*c_dag('do',0)*c_dag('up',2)*c('up',2)*c('do',0) + 0.8*c_dag('do',0)*c_dag('up',2)*c('up',0)*c('do',2) + 2.2*c_dag('do',1)*c_dag('do',2)*c('do',2)*c('do',1) + 0.8*c_dag('do',1)*c_dag('up',0)*c('up',1)*c('do',0) + 3*c_dag('do',1)*c_dag('up',0)*c('up',0)*c('do',1) + 0.8*c_dag('do',1)*c_dag('up',1)*c('up',2)*c('do',2) + 4.6*c_dag('do',1)*c_dag('up',1)*c('up',1)*c('do',1) + 0.8*c_dag('do',1)*c_dag('up',1)*c('up',0)*c('do',0) + 3*c_dag('do',1)*c_dag('up',2)*c('up',2)*c('do',1) + 0.8*c_dag('do',1)*c_dag('up',2)*c('up',1)*c('do',2) + 0.8*c_dag('do',2)*c_dag('up',0)*c('up',2)*c('do',0) + 3*c_dag('do',2)*c_dag('up',0)*c('up',0)*c('do',2) + 0.8*c_dag('do',2)*c_dag('up',1)*c('up',2)*c('do',1) + 3*c_dag('do',2)*c_dag('up',1)*c('up',1)*c('do',2) + 4.6*c_dag('do',2)*c_dag('up',2)*c('up',2)*c('do',2) + 0.8*c_dag('do',2)*c_dag('up',2)*c('up',1)*c('do',1) + 0.8*c_dag('do',2)*c_dag('up',2)*c('up',0)*c('do',0) + 2.2*c_dag('up',0)*c_dag('up',1)*c('up',1)*c('up',0) + 2.2*c_dag('up',0)*c_dag('up',2)*c('up',2)*c('up',0) + 2.2*c_dag('up',1)*c_dag('up',2)*c('up',2)*c('up',1)\n"
]
}
],
"source": [
"norb = 3\n",
"spin_names = ['up', 'do']\n",
"\n",
"U = 4.6\n",
"J = 0.8\n",
"\n",
"from triqs.operators.util.U_matrix import U_matrix_kanamori\n",
"\n",
"KanMat1, KanMat2 = U_matrix_kanamori(norb, U, J)\n",
"\n",
"from triqs.operators.util.hamiltonians import h_int_kanamori\n",
"\n",
"H_int_3 = h_int_kanamori(spin_names, norb, KanMat1, KanMat2, J, off_diag=True)\n",
"\n",
"from itertools import product\n",
"N_tot_3 = sum([ n(spin, idx) for spin, idx in product(spin_names, range(norb)) ])\n",
"\n",
"mu = 0.5*(5*U - 10*J)\n",
"\n",
"H_loc_3 = H_int_3 - mu * N_tot_3\n",
"\n",
"print(H_loc_3)"
]
},
{
"cell_type": "markdown",
"id": "aa879df5-ca9d-45fa-ae03-988074b80939",
"metadata": {},
"source": [
"## Step 2 - Spawn solver instance\n",
"\n",
"The second step in using `xca` is to instantiate an instance of the `triqs_xca.BlockSparseSolver` class."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "a210740f-d7db-4017-a81b-7ccd915f8758",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Starting run with 1 MPI rank(s) at : 2026-05-19 12:22:22.535403\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"____ ____________ _____\n",
"\\ \\/ /\\_ ___ \\ / _ \\\n",
" \\ / / \\ \\/ / /_\\ \\\n",
" / \\ \\ \\____/ | \\\n",
"/___/\\ \\ \\______ /\\____|__ /\n",
" \\_/ \\/ \\/ [github.com/TRIQS/xca]\n",
"\n",
"beta = 5.0, w_max = 50.0, eps = 1e-08, N_DLR = 32\n",
"AtomDiagReal: dim 4 with 4 subspaces dims [1] freq [4] E min/max +0.00E+00/+1.00E+00\n"
]
}
],
"source": [
"from triqs_xca import BlockSparseSolver as Solver\n",
"\n",
"norb = 1\n",
"\n",
"S = Solver(\n",
" H_loc=H_loc,\n",
" beta=5.0, \n",
" w_max=50.0,\n",
" eps=1e-8, \n",
" gf_struct=[('up', norb), ('do', norb)], \n",
" )"
]
},
{
"cell_type": "markdown",
"id": "8b4605f0-84d6-493b-a012-688506faf2c1",
"metadata": {},
"source": [
"The solver constructor takes the input\n",
"\n",
"- `H_loc` : Impurity local Hamiltonian\n",
"- `beta`: Inverse temperature\n",
"- `w_max`: DLR frequency cut-off (the spectrum of the model must be in the range `[-w_max, +w_max]`\n",
"- `eps`: Accuracy of Discrete Lehmann Representation (DLR) used for imaginary time response functions\n",
"- `gf_struct`: Green's function structure 1st index: name, 2nd index: dimension of subspace"
]
},
{
"cell_type": "markdown",
"id": "21baf8a1-51f3-4f38-b395-5eabd7ba1678",
"metadata": {},
"source": [
"### Using atomic symmetries\n",
"\n",
"The solver takes the optional parameter `conserved_operators` which defaults to `'automatic'`. This will automatically use the full set of state permutation symmetries that block diagonalizes the local Hamiltonian `H_loc`.\n",
"\n",
"**Note:** for multiorbital systems it can be faster to use only a subset of symmetries, like total density, by setting `conserved_operators=[N_tot]`, where `N_tot` is the total density operator (using Triqs 2nd quantization operators), i.e. `N_tot = sum([n('up', oidx) + n('do', oidx) for oidx in range(norb)])`. (To disable the use of symmetries simply pass `conserved_operators=[]`.)"
]
},
{
"cell_type": "markdown",
"id": "378098ff-7233-4434-9bec-9b75a31c6342",
"metadata": {},
"source": [
"## Step 3 - Hybridization function\n",
"\n",
"To fully specify the Anderson impurity model we also have to supply the hybridization function $\\Delta(\\tau)$ that describes the retarded fluctuation of electrons to the environment. The solver instance `S` sets up the hybridization function container available as `S.Delta_tau` and we need to explicitly set its value."
]
},
{
"cell_type": "markdown",
"id": "ff8a13fd-5eeb-46d1-9d5d-bf30825cf4af",
"metadata": {},
"source": [
"### Examples\n",
"\n",
"Here is a simple example with a Hybridization function $\\Delta(\\tau)$ comprised of a **single discrete pole/state**."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "0ada854e-a4b4-491c-afbd-5540c6d91a25",
"metadata": {},
"outputs": [],
"source": [
"from triqs.gfs import make_gf_dlr_imtime, make_gf_dlr_imfreq\n",
"from triqs.gfs import inverse, iOmega_n\n",
"\n",
"for bidx, delta_tau in S.Delta_tau:\n",
" delta_w = make_gf_dlr_imfreq(delta_tau)\n",
" delta_w << inverse(iOmega_n - 1.0)\n",
" delta_tau[:] = make_gf_dlr_imtime(delta_w)"
]
},
{
"cell_type": "markdown",
"id": "7ecac644-9fb2-43c0-b53d-8b26f055db75",
"metadata": {},
"source": [
"Another common test case is a spin- and orbital-diagonal hybridzation function with a **semi-circular density of states**."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "d2e4b357-677c-4238-a60a-dc3497f16648",
"metadata": {},
"outputs": [],
"source": [
"from triqs.gfs import make_gf_dlr_imtime, make_gf_dlr_imfreq, SemiCircular\n",
"\n",
"for bidx, delta_tau in S.Delta_tau:\n",
" delta_w = make_gf_dlr_imfreq(delta_tau)\n",
" delta_w << SemiCircular(2.0)\n",
" delta_tau[:] = make_gf_dlr_imtime(delta_w)"
]
},
{
"cell_type": "markdown",
"id": "27ec122a-487d-490f-9762-16236168f559",
"metadata": {},
"source": [
"## Step 4 - Self-consistent solution\n",
"\n",
"Combining the hybridization function `S.Delta_tau` and the local many-body Hamiltonian `H_loc` the Anderson impurity model is fully specified and we can call the `S.solve(...)` method to obtain the pseudo-particle self-consistent solution using all self-energy diagrams up to and including the expansion order `max_order`."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "e7c57247-8d19-4244-a16e-d9096bc3ea82",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"AAA: Error 1.87E-01 using 1 support and 30 fitting points (step 1/31)\n",
"AAA: Error 6.16E-04 using 2 support and 28 fitting points (step 2/31)\n",
"AAA: Error 2.40E-07 using 3 support and 26 fitting points (step 3/31)\n",
"AAA: Converged after 3 steps with error 2.40E-07.\n",
"TDC: Error 7.49E-07 for 3 AAA steps (Error 6.16E-04 no opt) c.f. tol 1.00E-05.\n",
"AAA: Error 1.87E-01 using 1 support and 30 fitting points (step 1/1)\n",
"TDC: Error 2.70E-01 for 1 AAA steps (Error 7.34E-01 no opt) c.f. tol 1.00E-05.\n",
"AAA: Error 1.87E-01 using 1 support and 30 fitting points (step 1/2)\n",
"AAA: Error 6.16E-04 using 2 support and 28 fitting points (step 2/2)\n",
"TDC: Error 1.58E-03 for 2 AAA steps (Error 1.87E-01 no opt) c.f. tol 1.00E-05.\n",
"TDC: Compression finished with 3 AAA steps and error 7.49E-07.\n",
"Adapol: Fit error = 7.49E-07 < tol_adapol = 1.00E-05, N_poles = 5\n",
"iter = 1, diff_G = 1.79E-01, Z-1 = -1.65E-10, eta = 5.14E-01\n",
"iter = 2, diff_G = 1.25E-01, Z-1 = -9.58E-11, eta = 6.81E-01\n",
"iter = 3, diff_G = 3.70E-02, Z-1 = -2.49E-11, eta = 7.23E-01\n",
"iter = 4, diff_G = 7.26E-03, Z-1 = -4.60E-12, eta = 7.31E-01\n",
"iter = 5, diff_G = 9.59E-04, Z-1 = -5.60E-13, eta = 7.32E-01\n",
"iter = 6, diff_G = 8.61E-05, Z-1 = -4.75E-14, eta = 7.32E-01\n",
"Converged after 6 iterations with diff_G = 8.61E-05 < tol = 1.00E-04\n",
"\n",
"Timing: incl. excl.\n",
"----------------------------------------------------\n",
"Adapol hybridization fit: 0.008 0.008 3.3% ||\n",
"AtomDiag Init: 0.000 0.000 0.0% |\n",
"DiagramEvaluator Init: 0.004 0.004 1.5% ||\n",
"Dyson: 0.054 0.001 0.4% |\n",
" Setup: 0.023 0.023 8.9% |---|\n",
" Solve: 0.030 0.030 11.5% |----|\n",
"Sigma: 0.135 0.000 0.1% |\n",
" Order 1: 0.004 0.004 1.6% ||\n",
" Order 2: 0.130 0.130 50.4% |-------------------|\n",
"Single-particle Gf: 0.012 0.000 0.0% |\n",
" Order 1: 0.000 0.000 0.1% |\n",
" Order 2: 0.011 0.011 4.4% |-|\n",
"Other: 0.046 0.046 17.8% |------|\n",
"----------------------------------------------------\n",
"Total: 0.259 100.0%\n",
"\n"
]
}
],
"source": [
"S.solve(max_order=2, tol=1e-4, verbose=True)"
]
},
{
"cell_type": "markdown",
"id": "3ce04ad5-b7a8-47a6-b52a-308659acfa7f",
"metadata": {},
"source": [
"### Hybridization compression\n",
"\n",
"For `max_order > 1` the solver will by default try to compress the hybridization function using [Adapol](https://github.com/flatironinstitute/adapol) to a minimal set of poles. The fit tolerance `hyb_tol` defaults to `0.1*tol` but can be passed as an optional argument to the `solve(...)` function.\n",
"\n",
"It is also possible to turn off the compression entirely by passing `hyb_comp=False` which will cause the solver to use the (larger) DLR pole representation, i.e."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "c367cfd4-61e6-477c-b9df-25d8a0b6a15f",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Hybridization: using DLR expansion with N_poles = 32\n",
"iter = 1, diff_G = 1.06E-04, Z-1 = -2.53E-14, eta = 7.32E-01\n",
"iter = 2, diff_G = 1.03E-05, Z-1 = +1.40E-14, eta = 7.32E-01\n",
"Converged after 2 iterations with diff_G = 1.03E-05 < tol = 1.00E-04\n",
"\n",
"Timing: incl. excl.\n",
"----------------------------------------------------\n",
"Adapol hybridization fit: 0.018 0.018 3.1% ||\n",
"AtomDiag Init: 0.000 0.000 0.0% |\n",
"DiagramEvaluator Init: 0.004 0.004 0.7% |\n",
"Dyson: 0.059 0.001 0.2% |\n",
" Setup: 0.023 0.023 3.9% |-|\n",
" Solve: 0.034 0.034 5.9% |-|\n",
"Sigma: 0.346 0.000 0.1% |\n",
" Order 1: 0.005 0.005 0.9% |\n",
" Order 2: 0.340 0.340 58.5% |----------------------|\n",
"Single-particle Gf: 0.083 0.000 0.0% |\n",
" Order 1: 0.001 0.001 0.1% |\n",
" Order 2: 0.083 0.083 14.2% |-----|\n",
"Other: 0.071 0.071 12.3% |----|\n",
"----------------------------------------------------\n",
"Total: 0.581 100.0%\n",
"\n"
]
}
],
"source": [
"S.solve(max_order=2, tol=1e-4, verbose=True, hyb_comp=False)"
]
},
{
"cell_type": "markdown",
"id": "bba188a4-d19d-477f-bdf0-b135f5c5d01a",
"metadata": {},
"source": [
"## Step 5 - Single particle response function\n",
"\n",
"After convergence the solver computes the diagrammatic series for the single particle Green's function that is available as the member `S.G_tau` of the solver class instance `S`."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "eeb420d4-fc66-436b-8f4a-483f35715ded",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Green Function G composed of 2 blocks: \n",
" Green's Function G_up with mesh DLR imaginary time mesh of size 32 with beta = 5, statistics = Fermion, w_max = 50, eps = 1e-08, symmetrized = false and target_shape (1, 1): \n",
" \n",
" Green's Function G_do with mesh DLR imaginary time mesh of size 32 with beta = 5, statistics = Fermion, w_max = 50, eps = 1e-08, symmetrized = false and target_shape (1, 1): \n",
" \n",
"\n"
]
}
],
"source": [
"print(S.G_tau)"
]
},
{
"cell_type": "markdown",
"id": "83a9a807-d760-424b-a80f-1cb2c66ff1d0",
"metadata": {},
"source": [
"## Step 6 - Store solver to disk\n",
"\n",
"The solver is hdf5 serializable and can be written to disk using the `h5` module."
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "76b80b13-05a5-49c1-ae86-ed868d5c693a",
"metadata": {},
"outputs": [],
"source": [
"from h5 import HDFArchive\n",
"\n",
"filename = f'data_xca_result.h5'\n",
"with HDFArchive(filename, 'w') as A: \n",
" A['S'] = S"
]
},
{
"cell_type": "markdown",
"id": "b7c6ef11-3bcc-43d5-892e-f08389329da0",
"metadata": {},
"source": [
"## Step 7 - Read result from disk\n",
"\n",
"For later visualization we can now read the solver and its data from disk."
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "477ab542-bd45-4aa8-b3d2-1851e7d5ce3d",
"metadata": {},
"outputs": [],
"source": [
"with HDFArchive(filename, 'r') as A:\n",
" S = A['S']"
]
},
{
"cell_type": "markdown",
"id": "c493ba7c-f201-467c-aa05-87cbe9b23db4",
"metadata": {},
"source": [
"## Step 7 - Postprocessing and visualization \n",
"\n",
"The single particle Green's function can now be visualized using the standard Triqs plotting tools."
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "f202ea70-e583-4e8c-bb8d-6cf3564ad5fe",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"from triqs.gfs import make_gf_imtime\n",
"from triqs.plot.mpl_interface import plt, oplot, oplotr, oploti\n",
"\n",
"oplotr(S.G_tau['up'][0, 0], label=None)\n",
"oplotr(make_gf_imtime(S.G_tau['up'][0, 0], n_tau=400), label=None)\n",
"plt.ylabel(r'$G(\\tau)$');"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.12"
}
},
"nbformat": 4,
"nbformat_minor": 5
}